%% pde_stokes_t
clc;
clear;

index = 6;

syms x y t nu
switch index
    case 1
        u1 = (x.^2.*y.^2+exp(-y)).*cos(2*pi*t);
        u2 = (-(2/3)*x.*y.^3+2-pi*sin(pi*x)).*cos(2*pi*t);
        p = -(2-pi*sin(pi*x)).*cos(2*pi*y).*cos(2*pi*t);
    case 2
        u1 = pi*sin(2*pi*y).*sin(pi*x).^2.*sin(t);
        u2 = -pi*sin(2*pi*x).*sin(pi*y).^2.*sin(t);
        p = 10*cos(pi*x).*sin(pi*y).*sin(t);
    case 3
        u1 = 10.*x.^2.*y.*(2.*y-1).*(x-1).^2.*(y-1).*cos(t);
        u2 = -10.*x.*y.^2.*(2.*x-1).*(x-1).*(y-1).^2.*cos(t);
        p = (20.*x-10).*(2.*y-1).*cos(t);
    case 4
        u1 = 5*(x.^2-2*x.^3+x.^4).*(2*y-6*y.^2+4*y.^3).*cos(2*pi*t);
        u2 = -5*(y.^2-2*y.^3+y.^4).*(2*x-6*x.^2+4*x.^3).*cos(2*pi*t);
        p = 5*(2*x-1).*(2*y-1).*cos(2*pi*t);
    case 5
        u1 = (x.^2-2*x.^3+x.^4).*(2*y-6*y.^2+4*y.^3).*cos(t);
        u2 = -(2*x-6*x.^2+4*x.^3).*(y.^2-2*y.^3+y.^4).*cos(t);
        p = (x.^2+y.^2-1.5).*cos(t);
    case 6
        u1 = -sin(pi*x).^2.*sin(2*pi*y).*cos(t);
        u2 = sin(2*pi*x).*sin(pi*y).^2.*cos(t);
        p = (x.*y-0.25).*cos(t);
end

u1dt = diff(u1,t);
u1dx = diff(u1,x);
u1dy = diff(u1,y);
u1dxdx = diff(u1dx,x);
u1dydy = diff(u1dy,y);
u2dt = diff(u2,t);
u2dx = diff(u2,x);
u2dy = diff(u2,y);
u2dxdx = diff(u2dx,x);
u2dydy = diff(u2dy,y);
pdx = diff(p,x);
pdy = diff(p,y);

f1 = u1dt - nu*(u1dxdx+u1dydy) + pdx;
f2 = u2dt - nu*(u2dxdx+u2dydy) + pdy;

disp("===u1导数===");
matlabFunction(simplify(u1dx), "Vars", {x,y,t})
matlabFunction(simplify(u1dy), "Vars", {x,y,t})
disp("===u2导数===");
matlabFunction(simplify(u2dx), "Vars", {x,y,t})
matlabFunction(simplify(u2dy), "Vars", {x,y,t})
disp("===p导数===");
matlabFunction(simplify(pdx), "Vars", {x,y,t})
matlabFunction(simplify(pdy), "Vars", {x,y,t})
disp("===非定常stokes方程===");
matlabFunction(simplify(f1), "Vars", {x,y,t,nu})
matlabFunction(simplify(f2), "Vars", {x,y,t,nu})